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Abstract 

This paper presents the development of a mesh adaptation 
module for a multilevel Cartesian solver. While the module 
allows mesh refinement to be driven by a variety’ of different 
refinement parameters, a central feature in its design is the 
incorporation of a multilevel error estimator based upon 
direct estimates of the local truncation error using r-extrapo- 
lation. This error indicator exploits the fact that in regions of 
uniform Cartesian mesh, the spatial operator is exactly the 
same on the fine and coarse grids, and local truncation error 
estimates can be constructed by evaluating the residual on 
the coarse grid of the restricted solution from the fine grid. A 
new strategy for adaptive h-refinement is also developed to 
prevent errors in smooth regions of the flow from being 
masked by shocks and other discontinuous features. For cer- 
tain classes of error histograms, this strategy is optimal for 
achieving equidistribution of the refinement parameters on 
hierarchical meshes, and therefore ensures grid converged 
solutions will be achieved for appropriately chosen refine- 
ment parameters. The robustness and accuracy of the adapta- 
tion module is demonstrated using both simple model 
problems and complex three dimensional examples using 
meshes with from lCr to / if cells. 

1 Introduction 

While error estimation using local gradient recovery tech- 
niques have long been popular in structural mechanics and 
other disciplines governed by elliptic systems, 1 such rigor- 
ous error estimates are generally not available in fluid 
mechanics where the governing equations contain non-self- 
adjoint operators. As a consequence, error estimation and 
adaptation for fluid mechanics has evolved over a different 
path. The simplest methods in the literature are gradient, and 
undivided difference based feature detectors> 11 While 
these methods have been extremely successful for various 
classes of problems 71 17 their lack of formalism makes them 
difficult to blindly apply to problems far from established 
experience. Indicators based upon estimates of interpolation 
error, 171181191 lend some of the missing formalism. However, 
since these essentially compute the local curvature of a repre- 
sentative variable or combination of variables, they are not 
clearly superior (or different from) straight feature detec- 
tion.^ 1 The multilevel error estimators from the literature on 
Adaptive Mesh Refinement (AMR) were among the first 
Richardson extrapolation- like error estimators 1 !UJ1 but 

have been only narrowly used outside the AMR community. 
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More recently, there has been interest in multilevel error indi- 
cators which use the residual of a higher order intemolation 
of the existing solution to estimate the local error} 

This approach is attractive since it can be combined with 
solution of an adjoint problem to produce a mesh adaptation 
strategy that seeks to minimize error in a specific integral 
quantity of engineering interest/ 13 ^ The hope of these 
techniques is that, for the cost of the adjoint solution and 
Jacobian storage, one can produce a mesh that focuses on 
reducing only those errors which are important to the prob- 
lem at hand, rather than simply equidistributing the error. 

While the research on adjoint-based mesh optimization is 
very exciting, two drawbacks make it currently unattractive 
of a general purpose analysis code. First, if the Jacobians are 
not already present in the solver, it is an expensive proposi- 
tion to form them solely to drive the adaptation. Secondly 
and more fundamentally, many problems of interest have 
competing requirements that cannot always be encapsulated 
into a single output functional. For example, a user may wish 
to optimize a single simulation for lift, drag and pitching 
moment. While this is currently an area of research interest, 
results are not currently in-hand. 

In this paper, we develop a new multilevel adaptation module 
for use on adapted Cartesian meshes with embedded bound- 
aries. Figure 1 shows a sample hierarchy of such meshes used 
by a new, parallel multigrid solver/ 1 ^ the fact that these 
coarse meshes can be automatically generated makes multi- 
level error estimation a viable alternative, even for grids 
around very complex geometries. Seeking to take advantage 
of this fact, we explore Richardson extrapolation-like error 
estimators which are modified to exploit the hierarchical 
nature of adapted Cartesian grids. Even without the sensitiv- 
ity information provided by adjoint-based approaches, the 
module will still offer huge savings over a priori mesh 
enrichment, and taking advantage of this is an important step. 
The estimators we present are light-weight, and both fast to 
implement and execute and leverage much of the machinery 
already built into multigrid solvers. 

Since Cartesian grid refinement is restricted to cell subdivi- 
sion (/z-refinement), the approach is intrinsically discrete. 
There exists a necessary reliance upon a refinement threshold 
which selects cells destined for subdivision. This work exam- 
ines the topic of threshold selection in detail. We develop a 
methodology for robustly setting this threshold which also 
ensures consistency. The approach both controls the level of 
error in the domain and equidistributes the remainin error as 
fast as possible within the restrictions of hierarchical (nested) 
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Figure 1 : Hierarchical sequence of partitioned meshes 
flow solver of ref. [17]. 

refinement. Given a robust strategy for setting the adaptation 
threshold, adaptive cases may be run fully automated. This 
reduces the sensitivity of the results to user experience/inter- 
vention which has long been a drawback of adaptive meth- 
ods. 


around lifting body configuration used by parallel 

satisfied the difference equation, the exact solution will not, 
and the difference will be the local truncation error, 
LTE hJc (x,t). 

LTE h k (x, t) = J[k(x, t + k)- u(x, /)] + *<“(*> -£ 2 ) (3) 


2 Multilevel Error Estimation 

Multilevel Richardson-type error estimators like those in [12] 
are attractive because of their low cost and firm theoretical 
underpinnings in smooth solutions. After developing these 
for perfectly nesting Cartesian mesh hierarchies, this section 
goes on to examine their application at adaptation boundaries 
and cells cut by geometry. It goes on to present shock and 
vortex-core detection strategies also planned as options to 
directly adapt to under-resolved features. 

2.1 Local TVuncation Error Measurement 

Consider the 1-D wave equation u t + A u x = 0 discretized 
with forward Euler in time and arbitrary spatial differencing 
using a timestep k , and a cell dimension h. 

if: 1 is the discrete approximation to the continuous solution, 
u(x, 0, in cell j at time n, while F L] and F R j are the numeri- 
cal fluxes through the left and right boundaries of j and con- 
tain the details of the spatial operator. The difference of the 
numerical fluxes is the residual and the discrete equation can 
be written more compactly as 

where the discrete residual operator /?(•) now contains all 
details of the numerical flux balance for cellj. 

The local truncation error , LTE, measures how well the dis- 
crete equation models the actual PDE throughout the domain. 
It is defined by replacing the approximate solution U* with 
the exact solution u(x , t) in equation (2). Since U” exactly 


If the problem reaches a steady state, then 
u(x,t+k) = u{x,t) and we can drop the dependence on 
time. 

LTE t (x) = (4) 

Equation 4 is instructive since it relates the local truncation 
error to the discrete residual operator R(*) and the cell dimen- 
sion h . Furthermore is clear from the derivation that in multi- 
ple dimensions LTE(x) ~ Residual/ h d where d is the number 
of spatial dimensions of the domain. 

22 Local Truncation Error on Uniform Cartesian 
Meshes 

Equation 4 is useful since it permits direct measurement of 
the local truncation error of a particular numerical method on 
an actual computational grid for any problem that has an 
exact solution. One such example is the supersonic vortex 
model problem used in reference 18. Figure 2 outlines the 
problem and shows 4 sample telescoping meshes. The actual 
meshes used had from 130 to 7800 control volumes. After 
initialization with the exact solution, the local truncation 
error was computed within each cell by dividing the residual 
by the cell volume (following §2.1 and eq. (4)) using the 
Euler solver in reference 17 without limiters. 

The norm of the LTE in density was computed for each 
mesh and is displayed in figure 3 as a function of the number 
of cells in each mesh. Performing a regression analysis on the 
Cartesian mesh data reveals an asymptotic slope of 2.1 1 . For 
smooth solutions, one can show that the global error is of the 
same order 19 . 
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For comparison, figure 3 also contains the results of the same 
experiment performed with three types of body-fitted grids: 
(1) regular, nearly unit aspect ratio quads, (2) right triangles 
made by subdividing the quads, and (3) a “quality” triangula- 
tion in which no angle was less than 29°. All methods 
used the same numerical scheme. While (he details of this 
experiment are contained in the Appendix, figure 3 displays 
the results. The asymptotic slopes of all four mesh types are 
shown in the figure. We note that the slopes of all but the 
“quality” mesh are similar but that the level of error on the 
Cartesian grids is consistently 6 to 10 times lower than the 
other grids over the full range of problem sizes. 



Number of control volumes 


Figure 3: L x norm of LTE in density for supersonic vortex 
model problem on sequence of meshes for 4 mesh. The 
asymptotic slope of each curve is labeled. Further 
details of this experiment are presented in Appendix A. 


23 Estimation of LTE on Nested Multilevel Meshes 

Since the LTE can be computed for each cell in the domain, it 
provides an excellent basis for constructing an refinement 
parameter within each cell. Unfortunately, since few prob- 
lems of practical interest have an exact solution, we must 
develop a method of estimating the LTE within each cell. 

In §2.2, the fact that we could reasonably draw straight lines 
through the data in figure 3 established the local and global 
order of the method. We can express this more formally by 
recalling that for a discrete method with local order p, there 
exists a constant Ci te such that 

\LTE h {x)\±C lte hP. (5) 

Section 2.2 also established that for the solver in reference 
17 , p was around 2.1 1 when ||^| is the 1-norm. Furthermore, 
since the RHS of equation 5 vanishes for small h, the method 
is obviously also consistent. Consistency implies that when h 
is small, the discrete solution Uj ^ in cell) of a mesh will be a 
good approximation to the exact solution u(x ) , and this 
statement is the basis of our multilevel error estimation pro- 
cedure. 

Assume that the numerical method has converged on a fine 
grid with mesh spacing h so that the residual of the discrete 
solution will be zero for every cell j in the domain. 

R h (Uj) = 0 (6) 

Since the method is consistent, the current discrete solution is 
assumed to be close to the exact solution, U j h ^ u{x ) . An 
estimate of the local truncation error on a coarse grid with 
mesh spacing H can then be written by substituting the dis- 
crete solution on the fine grid for the exact solution in equa- 
tion 4. 


LTE/x) 



(7) 


Following equation 7 we obtain an estimate of the LTE on the 
coarse grid by first restricting the discrete solution on the fine 
grid Uj h to a coarser grid using the interpolation operator 
I h , and then evaluating the discrete residual of the restricted 
solution on the coarse grid. Since Uj ^ satisfied the numerical 
scheme on the fine grid, the restricted solution I^Uj h will 
not, in general, produce exactly zero residual on the coarse 
grid. Thus, in the same way as the exact solution provided a 
means of measuring LTE /, in §2.1 , the discrete solution on the 
fine grid provides a method of estimating the error on the 
coarse grid. Since adaptively refined Cartesian meshes nest 
exactly, the residual operator R is exactly the same on coarse 
and fine meshes*. This permits us to apply equation 5 on a 
cell-by-cell basis. The coarse grid estimates are then used to 
trigger cell refinement on fine grid. This method of error esti- 
mation is known as t -extrapolation within the multigrid com- 


1 . In ref. 12 an estimate similar to eq.7 was used on curvilin- 
ear meshes. However, on curvilinear meshes, R /, * R h 
since the coarse grid control volumes are not the exact 
union of the fine grid control volumes. 
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munity/ 21 ^ and the perfectly nested residual operators that 
occur on Cartesian grids makes its use extremely attractive 
here. 

In the present work, the restriction operator /f is a volume 
weighted average, which is the same as that used by the mul- 
tigrid smoother used for convergence acceleration. By using 
the prolongation operator from the multigrid scheme as well, 
we are able to construct refinement parameters on the fine 
mesh for the cost of one restriction, one coarse mesh residual 
evaluation and one prolongation. Since these operations use 
machinery which already exists in the multigrid solver, they 
may be implemented with little effort. 

2 A Mesh Irregularities 

The preceding section examined error estimation on hierar- 
chies of uniform Cartesian meshes. Adaptively refined Carte- 
sian meshes, however, introduce some real-world 
complications which need discussion. As mentioned in the 
final paragraphs of §2.3, the multilevel error estimation pro- 
cedure relies upon use of the same spatial operator R on both 
the fine and coarse meshes. We denote the spatial operator on 
the pure Cartesian regions of the mesh as R ± . While R ± is 
clearly in use in uniform regions of the grid, nearer to adapta- 
tion boundaries and cut-cells the situation merits closer 
examination. 

Figure 4 depicts a sample fine grid (a) and the corresponding 
coarse cells (b) near a simple refinement boundary. The adap- 
tation boundary introduces irregularities into the connectivity 
graph which locally changes the residual operator. Thus, on 
the fine mesh, (Fig. 4a) all the cells which are immediately 
adjacent to the adaptation boundary have R * R L . The situa- 
tion is the same on the coarse mesh (Fig. 4b) and since the 
form of R h and that of R H must be identical for equation (5) 
to hold, the multilevel error estimator from §2.3 can’t be 
accurately constructed for the fine cells shown shaded in fig- 
ure 4a. Since this difficulty is associated with mesh irregular- 
ity, it appears in the residual operator of cut-cells as well. 

While the formal basis for local truncation error estimation 
makes it preferable to straight feature detection, there are 
several drawbacks that make them difficult to use - even for 
Cartesian multilevel meshes. First, for a reliable error esti- 
mate, the same stencil needs to be applied to both the coarse 
and fine grids. Thus, at mesh irregularities such as grid inter- 
faces and cutcells, the procedure described above does not 
produce a reliable estimate. Colella el. al [22] present a nice 
treatment of this problem in which the estimate is multiplied 
by the number of cells of that type. In essence, this scales a 
cell’s contribution to the error by the number affected cells. 
While the local truncation error at an interface is larger than 
at a uniform cell of the finer or even the coarser level, related 
work on convergence theory for non-uniform grids suggests 
that the estimate is still overly pessimistic. p3][24] (These 
results point out that its possible for 0(h) LTE to still pre- 
serve 0(h) global accuracy.) At present, we simply ignore 
error estimates at interfaces and cut cells, and rely on a buff- 
ering procedure to supply error values using piecewise con- 



Figure 4: Illustration of difference stencils against adapta- 
tion boundary between cells of different levels. Cells in 
the shaded region of the fine mesh restrict their solution 
into the corresponding coarse mesh shown. Since the 
residual operator on this coarse mesh is different from 
that in the corresponding cells on the fine mesh, a cell- 
by-cell interpretation of eq.5 is not valid and accurate 

stant extrapolation of the error estimates from neighboring 
cells. 

Inspection of figure 4 demonstrates that this irregularity 
affects a larger portion of the mesh than one may initially 
suspect since the cell next to an interface has its gradient cal- 
culation affected. This cell in turn contaminates the stencil of 
a cell two away from the interface. Even a minor change like 
this can be important if it occurs in a region of near zero mea- 
sured LTE. Additionally, since it is the coarse grid which pro- 
vides the estimates, the pollution can extend as far as 4 cells 
from the interface on the fine grid, so that the fine grid cells 
have no reliable estimate of their own ( cf . fig. 4). To produce 
grids with a stencil of 5 regular cells in each direction, large 
regions of uniformly refined cells need to be generated in the 
mesh. Although x-estimates have been used with great suc- 
cess in block structured AMR methods, they are somewhat 
more difficult to implement on unstructured multilevel Carte- 
sian. Asymptotic arguments contend that irregular cells are 
confined to only 0(N 2 ) cells in a 3-D mesh with 0(N 3 ) cells. 
In practice however, the difference between 0(N 2 ) and 0(N ) 
is not as great on the coarse mesh as it is for the fine mesh, 
and even reasonably refined meshes can have as much as 
30% of the cells contaminated by irregularity. When the 
neighbors of cells directly affected are also counted, esti- 
mates on over half the mesh may have some contamination. 
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A second problem with LTE estimates is the need for a some- 
what reliable solution to have already been computed. The 
approximate solution should be in the asymptotic region for a 
valid error estimate. This observation leads to a somewhat 
different strategy than that followed by feature detection 
strategy. Feature detectors generally perform well even when 
starting from a very coarse grid, making possible many 
cycles of adaptive refinement. For the LTE based refinement, 
our strategy is to start with a fairly refined initial gnd, and 
only use two or three additional cycles of LTE driven adap- 
tive refinement. In the examples here the initial grid uses 
geometry based refinement only, but initial grids generated 
from feature-detection based refinement could also be used. 

25 Prespecified Adaptation Regions 

Both LTE estimates and feature detection approaches suffer 
from a difficulty in localizing them to regions of interest. For 
example, errors in the wake behind a wing persist for quite a 
distance downstream. However, if a user is interested in the 
pressure distribution on the wing, research suggests that there 
is no need to refine the wake region very much. In the 
present work we simply allow the user to define a prespeci- 
fied (Cartesian) region in which adaptation is permitted. 
While not particularly elegant, this approach is adequate in 
practice. 

2.6 Feature Detection 

Feature detection attempts to use the current discrete solution 
to determine where mesh needs enhancement. Common 
schemes use first or second-order undivided differences and 
gradient information of various flow quantities. They attempt 
to “smooth out” computational space in hopes of producing a 
uniform error distribution. 

Since the link between flow features and truncation-error is 
not as formal as the Richardson-like LTE estimates discussed 
earlier, approaches in the literature vary significantly and 
everything from gradients and second unsca ^ ec ^ 

and scaled differences have been used. [2]l 

In 1991 Warren et ai showed that since gradients and second 
derivatives stay approximately constant with mesh refine- 
ment, they make poor refinement parameters. A cell with a 
high gradient will continue to have a high gradient even after 
subdivision. From the standpoint of grid convergence, we 
note that if the refinement parameter is to act as a substitute 
for the real LTE it should have the same asymptotic behavior 
as the LTE in smooth regions of the flow. For a p th order 
scheme, this means that halving the mesh should reduce the 
error by 2 P . 

For a first difference based quantity to mimic the LTE behav- 
ior of a second-order scheme, it must be scaled by the local 
mesh size, h. In one dimension, this leads to a first difference 
based refinement parameter of the form: 


J J i J tyj 


( 8 ) 


which is a normalized version of the first difference parame- 
ter advocated by Warren et alP® 

Similarly, a second difference based parameter should remain 
un-divided to give the same behavior/ 


* 1 2 4>j _ ♦/+ 1/2 •*•*/- 1/2 _ h 2 YA (9) 

T j 4> y 4»j 3 4 h 

The 1-D refinement parameters in eqs.(8) and (9) can be 
extended using a finite volume approach. This produces a 
vector refinement parameter, with components for each of the 
cell’s dimensions. The first difference based refinement 
parameter in the k ,h direction of cell j is: 


r j,k = h j,k 


2 (?<{>,• * k) 


♦/ 


( 10 ) 


where k is the unit vector in the ^ direction. 


These vector refinement parameters can be used to drive both 
isotropic and anisotropic cell subdivision as discussed in 
ref. [4], and a similar approach may be used for the second 
difference based parameters. 

Popular choices of <|> are the density, velocity magnitude or 
sometimes the local static pressure. In addition, combinations 
of these scalars are also possible. The investigations in sec- 
tion 4 use either density or velocity magnitude for detection. 


3 An Optimal Strategy for ^-Refinement 

The adaptation strategy seeks to refine the mesh using the 
LTE estimates or other refinement parameters from §2 to 
improve the solution. Unlike unstructured or structured body- 
fitted approaches, /--refinement (redistribution) is not an 
option with Cartesian grids, and cell sub-division (//-refine- 
ment) is the only alternative. The adaptation strategy takes 
the refinement parameters from section 2 as input and returns 
a boolean tag for each cell in the fine mesh. Subdividing each 
tagged cell will improve the accuracy of the discrete solution 
by reducing the local errors. 

Algorithm A outlines the adaptation procedure. 

Algorithm A: Adaptation Strategy 

Input: Vector of refinement parameters for each 

cell on fine mesh, r 

Output: Vector of cells tagged for //-refinement t . 

A.l Normalize: Normalize the refinement parameter vec- 
tor. r-j = |ry|/||HI® 


1 . Mesh optimization methods like those in references [15] 

and [16] localize the adaptation automatically when the 

user constructs a functional that defines the single output of 
interest. 


2. Eq. (9) differs from that presented in ref.[26] (eq.15) since 
their parameter is re-scaled by the local mesh dimension 
and will therefore vanish faster than the LTE as the mesh 
size is decreased. 
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A 2 Tag: Apply the adaptation criteria a( m ) to the normal- 
ized refinement parameter to produce a vector of 
tagged cells, x - a(r) 

A3 Rules: Modify the set of tagged cells to ensure the 
validity and smoothness of the output mesh. 

A.3.1 Buffer: Add buffer layers of tagged cells. 
r +~r b {x) 

A32 Smoothness : Filter for island/void suppression, 
x — r/x) 

A. 3 .3 Validity: Ensure adaptation boundaries do not 
exceed 2:1. x «- r 7 (x) 

A. 4 Output final vector of tagged cells x for subdivision. 

The final set of tagged cells output in A. 4 is then /i-refined 
and the solution vector is initialized on the new fine mesh 
with the prolongation operator from the multigrid scheme. 

In the following sub- sections we present a new strategy for 
reducing the refinement parameters to some predetermined 
level. The method is optimal in that it accomplishes this task 
as fast as possible. In addition the method ensures that the 
permissible level of the refinement parameter chosen at the 
outset will not vary as the mesh and solution evolve. 


3.1 Equidistribution of Refinement Parameters 

Equidistribution aims at producing a mesh which contains 
the same level of refinement parameter in each cell. Since the 
refinement parameters are stand-ins for the local truncation 
error, this goal spreads the remaining discretization error out 
evenly over the domain. As this level is reduced, the method 
is guaranteed to converge to the correct solution. This princi- 
ple guards against over-resolving some features of the flow 
while overlooking others. 


In practice, equidistribution is somewhat over conservative 
most of the time. It assumes that all errors are equally impor- 
tant to the simulation, and this is certainly not the case most 
of the time. However, without additional guidance about 
what is important for a particular simulation, equidistribution 
simply ensures that everything in the simulation is equally 
correct. In addition, if a method can control the LTE distribu- 
tion to achieve equidistribution, then it can control the error 


Ncells-\ 



% 


Refinement Parameter 


Figure 5: Idealized histogram of refinement parameters 
in a mesh which as achieved equidistribution. 


to achieve a different goal. Its easy to conceive of inverse-dis- 
tance weightings or error weightings that take the local char- 
acteristics into account in order to identify those errors that 
have the strongest impact the output functionals of interest. 

Figure 5 shows the histogram of refinement parameters in a 
mesh which has achieved equidistribution, since all cells in 
the domain have the same error, they fall in the same bin, and 
the histogram looks like a delta function whose height is 

N cells' 

Like most strategies for adaptation, the paradigm of Alg. A is 
to start with some initial distribution of refinement parame- 
ters, and drive this distribution toward the idealized distribu- 
tion shown in fig. 5. 

Figure 6 shows the Gaussian-like distribution of refinement 
parameters which serves as a model for the histogram prior to 
/j-refinement. A common approach found in the literature is 
to set the refinement threshold to some fraction of a standard 
deviation above the mean of the distribution P**- 

32 A Fresh Look at Refinement Histograms 

Figure 7 shows an actual refinement histogram resulting from 
a coarse grid simulation (3775 cells) of flow over an ONER A 
M6 wing at transonic conditions^ This plot bears little 
resemblance to the idealized Gaussian-like model shown in 
figure 6. The refinement parameters lie between 0 and 1 . The 
mean value is 0.011 but the standard deviation 0.04. More- 
over, fully 82% of the cells lie below the mean, and almost 
50% have \r\ < 0.001 . As a consequence of this extreme dis- 
parity in scales, setting the threshold, f, anyplace above the 
mean only addresses the error in a handful of the most severe 
cells. Only after the very worst errors are reduced by many 
cell refinements will error in the bulk of the domain be 
addressed. In shocked flows, the refinement parameter will be 
highest in cells with shocks or other strong non-linear fea- 
tures. As a consequence, this approach will inadvertently 
result in a refinement process which over-resolves shocks and 
other severe features without ever addressing smooth regions 
of the flow. Oversights such as these have been shown to pro- 
duce an adaptive procedure which can actually converge to 
the wrong solution^. 



Figure 6: Idealized distribution of refinement parameters 
prior to ^-refinement. Region a contains N a cells. 
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Figure 7: Histogram of adaptation parameter for a 
coarse mesh simulation of flow around an ONERA 
M6 wing at transonic conditions^ 


In earlier work [4] we advocated a filtering approach which 
removed cells containing shocks or other strong features in 
an attempt to clean up the histogram prior to computing the 
mean and standard deviation. Even this approach, however is 
dubious given the huge disparity in scales. 

A less ad-hoc approach toward compressing the scales in the 
histogram is to simply take the log of the refinement parame- 
ter prior to binning up the histogram. Figure 8 shows the his- 
togram of the same data as in figure 7 but computed using 
log 2 (|r|) rather than simply |r| . The mean value of this new 
distribution is -6.4 and the standard deviation is 4.3. The 
rescaled data much more closely resembles the idealized dis- 
tribution in figure 5 . 



log_2(Rclincmcm Pantm.i 


Figure 8: Histogram of data from fig. 7 computed using 
log 2 (|r|) rather than the raw refinement parameter. 


33 Optimal Threshold Selection 

The motivation for choosing base 2 for the logarithms used to 
rescale the refinement parameters in the proceeding section 
becomes clear when selecting an appropriate threshold. Near 
grid convergence, each 2:1 cell refinement using a p -order 
scheme will reduce the LTE by a factor of 2 F . With 2:1 cell 
refinement, and the present second-order scheme, the chil- 
dren of a A-refined cell will therefore get translated an abso- 
lute distance of 2 units to the left on these base-2 histograms. 
Figure 9 illustrates this process. If there are N a cells in region 



■ Refinement Parameter 

Figure 9: /i-refinement moves the cells in region a to the left 
according to the order of accuracy of the scheme. If a 
contains N a cells, then a * contains m N a cells, where m 
is the number of children produced by refining a cell. 

a of figure 9, and each cell is subdivided into m children, then 
a* will contain m N a new cells. 

If our goal is equidistribution, then we desire to build the 
delta function of figure 5. An optimal method construct these 
as rapidly as possible. Assuming that the histogram is mono- 
tonically decreasing to the right of the median, the new histo- 
gram grows most rapidly if the highest point of a* is placed 
on top of the median of the existing distribution. Since the 
cells in a move 2 units to the left, the threshold which builds 
the highest peak is identically 2. 

Subsequent refinements will add to this same peak, and the 
target level of error will remain constant as the mesh evolves. 
It is interesting to note that if the threshold is chosen above or 
below than this value the target error level will continue 
migrate higher or lower (respectively) with subsequent cell 
refinements. 

Just as refinement moves cells to the left on the histogram, 
coarsening transfers cells to the right. In the absence of 
coarsening, these low-error cells will remain in the histogram 
and appear as a “tail” to the left of the peak value. Figure 10 
illustrates the evolving histogram. With the threshold chosen 
as described above, newly refined cells will not alter the his- 



Figure 10: Evolution of a histogram for a mesh without 
coarsening. Growth in the number of cells before the 
peak will be approximately geometric. 
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togram to the left of the peak value, and therefore no newly 
refined cells can ever end up in this tail. Since they were cells 
in the original unadapted mesh, this tail contains only coarse 
cells, and since coarse cells fill space very quickly, there can- 
not be very many of them in the domain as compared to the 
number of highly refined cells in the peak. In addition, since 
the peak was built via cell subdivision, and the number of 
children produced per cell is generally constant, the growth 
approaching the peak from the left must be approximately 
geometric. 

3.4 Specifying an Error Level 

Skeptics may point out that setting the adaptation threshold 
to its optimal value removes the user’s “control” of the adap- 
tive process. While this is precisely the goal of automation, 
there is a clear need for a user to be able to have some control 
over the level of error in the final solution. 

In three dimensions, isotropic /z-refinement of almost any 
polyhedron produces cells very quickly. Hexahedra are 
refined 8: 1 , and depending on the method of subdivision, tet- 
rahedra multiply at about the same rate. Even allowing for 
anisotropic subdivision of some cells, the 3-D growth rate is 
such that adaptively refined cells virtually always account for 
the median value in the histogram. 

Since the location of the peak can be controlled by adjusting 
t, the most efficient strategy is to set the desired peak location 
on as coarse a mesh as possible ( i.e . the first adaptation 
cycle). Subsequent adaptations will then continue to build on 
this peak using the optimal threshold. This allows the user to 
set a desirable error level based upon the histogram of the 
unadapted coarse mesh, and then drive the refinement hands- 
free. 





35 An Illustrative Example 

Using the coarse mesh simulation from the base-2 histogram 
in figure 8, figure 1 1 shows the evolution of the histogram in 
this simulation over the next 5 adaptation cycles. This exam- 
ple clearly shows the rapid growth of the peak in the histo- 
gram confirming its approach toward equidistribution of the 
refinement parameter. 

Since this is a real 3-D transonic flow, several issues merit 
discussion. Adaptation was driven by the undivided first dif- 
ference of velocity magnitude scaled by the local mesh 
dimension as presented in §2.6. The adaptive procedure in 
Alg. A tagged cells according to adaptation criteria, and these 
tags were then modified to satisfy the smoothness and mesh 
validity rules detailed in A3. In addition, the adaptation crite- 
ria used in this example was: 


a(rj) 


1 \rj\ > |r| + 2 

0 otherwise 


(ID 


Where \r\ is the mean value of the magnitude of the refine- 
ment parameter, rather than the precise median as called for 
by the theoretical development earlier in this section. Of 
course for the narrow peaked histograms shown in the figure, 




Figure 11: Evolution of histogram through 5 adaptation 
cycles for transonic ONERA M6 wing case using the 
optimal threshold. 
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the mean is close to the median - it is within one unit of the 
median at every cycle after the first. Nevertheless one would 
expect somewhat better performance if the true median was 
used. 

4 Numerical Examples and Discussion 

The LTE estimates and feature detection approaches in §2 
have been applied to both simple and complex configurations 
in three dimensions. This section begins by presenting results 
showing that, when combined with the /z-refinement strategy 
of §3, both can provide valid approaches to achieve grid-con- 
verged solutions. We then present adapted solutions on two 
complex configurations to demonstrate the robustness of the 
procedure when run “hands-off’ on real-world complex con- 
figurations. 

4.1 ONERA M6 Wing 

In ref. [17] the baseline Euler solver was validated using the 
well-known ONERA M6 wing example. 1251 This case con- 
siders the transonic flow over this wing at =0.84 and 
a = 3.06°. Although viscosity was obviously present in the 
experiment, the case has been widely studied using inviscid 
solvers, and a multitude of Euler solutions are available in the 
literature for comparison. 



Figure 12a: Flow over an ONERA M6 wing at M* =0.84 
and a = 3.06°, adaptation driven by x-estimates of 
LTE in density. The final mesh contains 1 .8M cells 


This flow was computed using both the multilevel x-extrapo- 
lation LTE estimates and scaled first-difference (eq.10) based 
feature detection using (j> = I n as a refinement parameter. 
Figure 12a displays both the mesh and solution resulting 
from the LTE based adaptation, while 12b contains these 
same plots using feature detection. In figure 12, the x-extrap- 
olation analysis used an initial mesh with 8 levels of geome- 
try based refinement and then an additional three cycles of 
solution-based refinement following the philosophy of §2.4. 
The final mesh shown contains 1 .8M cells. The feature detec- 
tion based results in figure 12b began on a mesh created with 
5 levels of geometry-based adaptation and about 30k cells 
followed by 6 cycles of adaptation. The final mesh contained 
1 ,9M cells. After 5 levels of adaptive refinement, the mesh 
contained about 1.2M cells and integrated quantities (normal, 
axial and side force) were virtually the same as those 
obtained on the final (1 .9M cell) mesh. The integrated quanti- 
ties for both examples changed by less than 0.1% in the last 
adaptation cycle suggesting that the results are grid con- 
verged. 

Lift and drag coefficients for x-extrapolation were: 0.3041 
and 0.01 17 while those the feature-detection were 0.3042 and 
0.0116. Comparison between the two simulations reveals a 
difference of less than 0.04% in the magnitude of the total 
force on the wing, and the final meshes have cell counts 
within 6% of each other. 



Figure 12b: Flow over an ONERA M6 wing at M x =0.84 
and a = 3.06°, adaptation driven using scaled first 
differences of velocity magnitude. The final mesh 
contains 1 .9M cells 
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Figure 13: Pressure profiles from transonic ONERA M6 
wing case at 44% span showing evolution of the C p 
history over the last 3 adaptation cycles. 

Figure 13 shows convergence of the C p profile at the 44% 
span station and includes an overplot of the experimental data 
for the simulation in figure 12b^ \ Behavior at this span sta- 
tion is typical and shows convergence of the adaptive process 
over the last 4 adaptation cycles. In the figure, the profiles of 
the solution after 5 and 6 adaptation cycles are essentially 
indistinguishable. Comparison with the experimental data is 
generally good and shows the same kinds of differences 
reported by other inviscid simulations of this viscous flow. In 
particular, the separation bubble following the rear shock is 
not modeled, and this shock is positioned slightly behind that 
in the experiment. 

42 Complex Configurations 

Figures 14 and 15 display the first of two examples showing 
real-world applications of the adaptation module on complex 
geometry. The supersonic canard-controlled missile geome- 



Figure 14: Geometry and adapted mesh for canard-con- 
trolled missile example. The final mesh has 4.5M cells 
and used 6 cycles of adaptive refinement, the last 3 of 
which were confined to the pre-specified adaptation 
box illustrated. 


try in figure 14a contains several features which make it chal- 
lenging to simulate. At angle -of- attack, or anytime the 
canards are deflected, they will create vortices which may 
interact with the tail fins of the missile. Due to the high fine- 
ness ratio of the missile, these vortices must convect over 30 
canard chord lengths before reaching the tail. Clearly excess 
numerical dissipation can easily destroy this important inter- 
action. In addition, the disparity in length scales on the geom- 
etry makes this simulation challenging. The canard chord is 
only 1 40 th of the body length. The simulation must resolve 

not only fine geometric scales like the leading and trailing 
edges of the canards, but also the bow shock on the missile, 
and the shock system generated by the canards themselves. 
Inviscid overset (structured) grid simulations with the Army’s 
OVERFLOW-D solver used over 30M points to resolve the 
features of this flow field. 

Supersonic flow over this missile was computed at zero 
degrees roll, and M 00 = 1,6 at a = 3°. The canards are 
deflected 15° (nose up). These conditions are give a reason- 
ably strong interaction between the canard tip vortices and 
the leeward pair of the interdigitated tail fins . 

Figure 16 shows contours of velocity magnitude in the dis- 
crete solution of this flow on the adapted mesh shown in fig- 
ure 14. The refinement parameter in equation (10) based on 
density was used to drive the adaptive process. Both figures 
clearly display the trajectory of the canard vortex system as 
they convect down the body. The final mesh has 4.5M cells 
and used 6 cycles of adaptive refinement, the last 3 of which 
were confined to the pre-specified adaptation box illustrated. 
Several axial cutting planes in figure 15 display the evolution 
of the canard vortex as it travels down the missile body. The 
computed normal force coefficient on the final mesh matches 
the inviscid results in ref [27] to within 3%. 

One interesting aspect of the missile simulation is that with 
the adaptive strategy outlined in §3 and the refinement 



Figure 15: Velocity magnitude contours for flow over a 
canard-controlled missile of in fig. 14 at 1.6 at 

a = 3° with the canards deflected 15° (nose up. 
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parameter from equation (10), the canard vortex and other 
important, but smooth regions of the flow are refined to the 
same level as the shocks. Thus the case can be made that the 
shocks are not receiving undue attention from the refinement 
scheme. 

This observation is further supported by the space-shuttle 
configuration displayed in figures 16-18. In this case the 
model was composed of 22 separate components, and 
included spoilers, flaps, rudders, engine bells, and other geo- 
metric detail .While the elevons and spoilers are nominally 
undefeated, some gaps exist, and there is flow leakage near 
these control surfaces. 

The half-body, power-off configuration was simulated at 
Moo= 1.5 and a = 8°, and refinement was focused in a box 
extending 3 body lengths in the crossfoot directions truncated 
just downstream of the orbiter. Figures 16-18 all display the 
solution using contours of velocity magnitude, and the 
scaled, undivided first difference of density was used as the 
adaptation parameter. Five cycles adaptation were earned out 
(hands-off) from an initially geometry refined mesh, produc- 
ing a final mesh with 8.5M cells. Figure 16 shows a nose-on 
view of the grid mirrored to the starboard side, with contours 
of the solution on the port side. While the bow and wing 
shocks are reasonably well resolved, nearly equal mesh spac- 
ing is used in much of the near-body lee-side flow. The flow 
structure is reasonably complex, with the body, wing, OMS 
pods, etc. all producing massive curved shocks that are cap- 
tured by the refinement. The curvature of these shocks results 
in strongly non-linear flow downstream and the refinement 
extends into these regions. In addition, both the wing tip vor- 
tex, and a vortex emanating from the gap between control 
surfaces are evident in this figure. 



Figure 16: Computational mesh and velocity contours of 
solution for orbiter simulation at M& = 1 .5 and a = 8°. 
The final mesh contains 8.5M cells. 


Figure 17 provides additional insight, displaying both the 
mesh and solution from a vantage point behind and above the 
port wing. The cutting plane in this view is located mid-way 
over the starboard wing, and some gaps between the control 
surfaces are visible. The bow, canopy, wing and trailing edge 
shocks are all clearly visible in this view. Figure 1 8 is a pro- 
file shot which contains a cutting plane through the symmetry 
plane to provides a better view of the canopy and bow 
shocks. 



Figure 17: Rear three-quarter view of orbiter geometry and 
mesh showing gaps between control surfaces and cut- 
ting plane through solution at a mid- span location. 
1.5, a = 8°, velocity contours 
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Figure 18: Side view of orbiter simulation and symmetry 
plane solution. M ac = 1.5, a = 8°, velocity contours. 


5 Conclusions and Future Work 

We have presented Cartesian mesh adaptation strategies 
driven by either local truncation error estimates or feature 
detection. The adaptation module adds a solution-based mesh 
refinement capability to the geometry based refinement of the 
Cartesian mesh generator. Timings have shown the new algo- 
rithm to be a bit faster than the original mesh generator, due 
to a few improvements in the algorithms and implementation. 
Both simple studies and highly complex 3D examples were 
presented with very high resolution, demonstrating the 
robustness and utility of the adaptation. The module produces 
several million Cartesian cells-per-minute on desktop com- 
puters and was demonstrated on complex example geome- 
tries with ~10 7 cells. 

An interesting highlight of this work is an optimal strategy 
for h-refinement based on log 2 ( ) histograms. This strategy 
avoids many of the pitfalls of the mean and standard devia- 
tion based approaches found in the literature. For hierarchal 
meshes, the approach is optimal in that it maximally equidis- 
tributes the refinement parameter for a given number of adap- 
tation cycles. We believe it to be more reliable and robust 
than mean and standard deviation based approaches. 

Our initial investigation of multilevel local truncation esti- 
mates was disappointing due to the irregularities in the mesh 
at the embedded boundaries and interfaces. While these com- 
plications can be overcome, they affect more of the mesh 
than was initially expected. More investigation of this prob- 
lem is needed. When applied to a model problem with a 
known analytic solution, these truncation error estimates re- 


confirmed the accuracy advantages in both LTE and global 
error offered by Cartesian meshes. 

Future work will focus on several outstanding topics. The 
behavior of the refinement strategy needs to be validated over 
a wider range of input conditions. While this strategy has 
been performed extremely well in initial investigations of 
flows with free stream Mach numbers from 0.8-1 .6, it has not 
been exhaustively tested for broader conditions. The strategy 
assumes that the histogram of refinement parameters is 
monotonically decreasing to the right of the median bin, and 
the validity of this assumption should be investigated under 
more extreme conditions. A second topic for further investi- 
gation is selection of an initial threshold to control the overall 
error level in the simulation. Currently this is done by inspec- 
tion of the coarse mesh histogram, but a more automated pro- 
cess would be desirable. We have also not implemented a 
mesh coarsening algorithm, and while this is not a major con- 
cern for steady flows, one will be needed for unsteady appli- 
cation in the future. Finally, adaptation is currently “focused” 
through the use of pre-specified adaptation regions. Such 
regions could be automatically generated using characteristic 
information from the flow to appropriately weight or 
unweight refinement parameters depending upon the cell’s 
location in the domain and the input Mach number. 
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8 Appendix 

The supersonic- vortex model problem discussed in §2.1 and 
figures 2 and 3 was performed to measure the local truncation 
error of the embedded-boundary Cartesian solver. In this 
appendix, we present details of this investigation and com- 
pare the Cartesian results with 3 other meshing systems. In 
addition to the nested Cartesian grids, this experiment was 
performed using three popular body- fitted meshing schemes, 
and the same underlying solver (upwind, with linear-recon- 
struction). Regular (structured) quad, right triangular, and 
“quality” triangular^ 81 body-fit meshes were used. Figure 
A.l shows the second-coarsest mesh of each of these types of 
grid. Four meshes of each type were used in the investigation 
and the meshes contained from 128 to 7809 control volumes. 
Special care was taken to match the numbers of control vol- 
umes for each mesh type as closely as possible. 

The example was computed with an inner Mach number, 
M in = 3.0, and taking p in = 1/y, p,„ = 1 , r, = 1 and r 0 = 1.9. 
Each case was initialized with the exact solution and the LTE 
(of density) within each cell was computed using eq.(4) by 
applying the residual operator without using flux limiters. 

Table A.l contains the L; norm of the LTE for each of the 16 
simulations. The simulations for each of the mesh types cor- 
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Triangles 


Figure A.l: Representative “Quality” Triangular, Right Tri- 
angular, and Regular Quad body- fit meshes used in LTE 
investigation. These meshes are the second coarsest 
used and have 505, 525 and 525 control volumes 
(respectively). 

related closely to a straight line, and the asymptotic slopes of 
each are given in the table. 


Some aspects of these data merit discussion. As noted in 
§2.2, the rate of convergence of the LTE for all mesh types 
are similar except for that of the “quality” triangular mesh. 
While all the other mesh types demonstrated second-order 
accuracy, results for this grid system were only slightly better 


than first-order. Since this mesh is not quite uniform, the 
stencils used for the gradient estimation and reconstruction 
on all of the “quality” meshes are all slightly irregular. Since 
every stencil is irregular, each is polluted to some degree by 
stretching, cancellation of errors cannot occur, and the result 
is a marked degradation in accuracy. 

These results contrast with earlier results for a similar prob- 
lem using nearly-equilateral triangles^ 18 ^. That investigation 
showed that regular equilateral triangles performed as well 
(or better) than regular quads or right triangles. In this case, 
however, the “quality” mesh is not equilateral, although all 
the triangles are well formed as guaranteed by the Ruppert’s 
2-D delaunay technique in ref. [28]. Quality 2-D meshes were 
chosen for this investigation since it is not possible to gener- 
ate uniform meshes of equilateral tetrahedra in 3-D. As a 
result, tetrahedral mesh generators typically resort to produc- 
ing “quality” meshes that guarantee some angle criterion is 
everywhere met, just as Ruppert’s Delaunay algorithm does 
in 2-D. 

The structured quad and right-triangular meshes are substan- 
tially smoother than the quality triangular meshes. Neverthe- 
less the LTE measurements indicate that even the mild 
irregularity in their stencil degrades their performance. While 
both provide second-order accuracy, the absolute error level 
is from 6 to 10 times worse than the Cartesian grid’s perfor- 
mance where irregularity is confined to the boundary. In fig- 
ure 3 the structured quads require nearly 10 times more cells 
to match the absolute level of LTE in the Cartesian scheme, 
and while the right triangular mesh performs slightly better, it 
does so at the cost of 50% more flux evaluations due to the 
higher edge count on these meshes. 


Cartesian Mesh with Embedded Boundary 

# of Control volumes Measured Lj { density ) Error 

138 0.03065 

507 0.00930 

1928 0.00246 

7549 0.00059 Asymptotic slope = 2.11 

Body-Fit Structured (Quad) Mesh 

# of Control volumes Measured LI { density ) Error 

144 0.30998 

525 0.09223 

2001 0.02422 

7809 0.00629 Asymptotic slope - 1.94 


Body-Fit Right Triangular Mesh 

# of Control volumes Measured L } (density) Error 

144 0.37926 

525 0.07571 

2001 0.01565 

7809 0.00347 Asymptotic slope = 2.28 

Body-Fit Quality Triangular Mesh 

# of Control volumes Measured L } ( density ) Error 

128 0.52552 

505 0.22529 

1918 0.11936 

7490 0.05940 Asymptotic slope = 1 .02 


Table A.l : Lj-Norm of LTE in density for each of the 16 meshes used in the supersonic vortex investigation. The 
“Quality” triangulation meshes were produced using the quality Delaunay triangulation algorithm of ref.[28] and 
had no angle less than 29° . 
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